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The article delivers the only still unknown coefficient in the 4th post-Newtonian energy expression 
for binary point masses on circular orbits as function of orbital angular frequency. Apart from a 
single coefficient, which is known solely numerically, all the coefficients are given as exact numbers. 
The shown Hamiltonian is presented in the center-of-mass frame and out of its 57 coefficients 
51 are given fully explicitly. Those coefficients are six coefficients more than previously achieved 
[Jaranowski/Schafer, Phys. Rev. D 86, 061503(R) (2012)]. The local divergences in the point-mass 
model are uniquely controlled by the method of dimensional regularization. As application, the last 
stable circular orbit is determined as function of the symmetric-mass-ratio parameter. 
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I. INTRODUCTION 

The measurement of gravitational waves will finally re- 
veal the motion of binary black holes with highest preci- 
sion including the 4th post- Newtonian (4PN) approxima- 
tion under consideration in the present article [l| . Partic- 
ularly important is the knowledge of the last stable circu- 
lar orbit (LSCO) beyond which binary systems do plunge 
and undergo merging. Though PN schemes are not quite 
appropriate for describing plunge orbits reasonable esti- 
mates can be made for the LSCO. For achieving more, 
the effective-one-body (EOB) approach by Damour and 
collaborators Q , which is strongly based on PN approxi- 
mations, allows treatment of the dynamics beyond LSCO 
(see, e.g., Of most urgent need is the knowledge of 

the PN dynamics beyond 3.5PN order, particularly the 
4PN order. 

In this communication we will report on the calculation 
of six more coefficients of the 4PN center-of-mass Hamil- 
tonian for binary point-mass systems which resulted in 
the knowledge of 51 coefficients out of total 57 ones. 
Without knowledge of the missing six coefficients, the 
full energy expression of binary systems on circular orbits 
is derived for the first time taking into account a coeffi- 
cient which was recently calculated numerically [f| . The 
numerically-only-known coefficient is uniquely connected 
with the global structure of the near zone and needs for 
its calculation the regularization of infrared divergences. 
All the other coefficients in the energy expression for cir- 
cular orbits at 4PN order - apart from the intimately 
connected with logarithmic function term linear in the 
ratio of reduced mass to total mass - are calculable with 
the aid of the technic of dimensional regularization, in 
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this way controlling the ultraviolet divergences resulting 
from the applied point-mass model. It is worth noticing 
that the corresponding logarithmic terms in the Hamilto- 
nian presented below were uniquely calculated with the 
aid of analytic regularization in three-dimensional space; 
here, pure dimensional regularization would have sup- 
plied logarithmic expressions not in agreement with the 
Poincare algebra. The previous achievements in the ex- 
pression for the energy on circular orbits at 4PN Q are 
confirmed through UH and, most recently, by [9j]. 

We employ the following notation: x = (V 1 ) (i = 
1,2,3) denotes a point in the 3-dimensional Euclidean 
space K 3 endowed with a standard Euclidean metric and 
a scalar product (denoted by a dot). Letters a and b 
(a, 6 = 1,2) are body labels, so x a G M. 3 denotes the po- 
sition of the ath point mass. We also define r a = x — x a , 
r a = \r a \, n a = r a /r a ; and for n / I), r ab = x a - x 5 , 
Tab = I r a b I, n ab = r ab /r ab ; \ ■ \ stands here for the Eu- 
clidean length of a vector. The linear momentum vector 
of the ath body is denoted by p a — (p a i), and m a de- 
notes its mass parameter. We abbreviate S (x — x a ) by 
<5 a . Extensive use has been made of the computer- algebra 
system Mathematica. 



II. THE ADM CANONICAL APPROACH IN d 
SPACE DIMENSIONS 



Let D = d + 1 denote the (analytically continued) 
space-time dimension. The ADM approach [l(| uses a 
d + 1 split of the coupled gravity-matter dynamics and 
works with the canonical pairs (x l a ,p a i) and ( gij , tt ij ) 
(i,j,k,... denote spatial indices taking formally d val- 
ues; a,b — 1,2 labels the particles). The dimensionally 
continued hamiltonian and momentum constraints read 
(in units where 167tGd = c = 1 with Gd denoting the 
generalized Newtonian gravitational constant and c the 
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speed of light) 



\9ik 9ji n ij ir M - -^j ( 9ij 



+ y2( m l+9a Par Paj )* 5 a , 



(2.1a) 



(2.1b) 



Here R denotes the space curvature of the hypersurface 
t = const, g is the determinant of the matrix (gij), 
S a denots <i-dimensional Dirac delta distribution (with 
/ d d x5 a = 1), g % i = 3 ij (Xu) is the finite part of met- 
ric evaluated at the particle position (which can be per- 
turbatively unambiguously defined), and Dj is the d- 
dimensional covariant derivative (acting on a tensor den- 
sity of weight one) . 

Taking into account the appropriate coordinate condi- 
tions 



where 



* = 1 



Id- 2 



Ad- 1 



(2.2) 



(2.3) 



holds and where hj^ is a symmetric transverse-traceless 
(TT) quantity, hj^ 1 = hj^ = 0, and the field momentum 
■k 1 ' 3 is splitted into its longitudinal and TT parts, respec- 
tively n ij = n ij +tt^ t , with 7T y = drf +dj7r l - \^d h -K k , 
from the constraint equations, the reduced Hamiltonian 
results in the form, 

fired [x a , p a , hj^ , 7T^ T ] = - J d d x A(j) [x a , p a , hj? , 7T^ T ] . 

This Hamiltonian describes the evolution of the matter 
and independent gravitational field variables. 

An autonomous conservative Hamiltonian can be ob- 
tained through the transition to a Routhian description, 

R [x a , p a , hjy, hj?] = H Ied - J d d x /'7; 1 - (2.5) 

Then the matter Hamiltonian reads 

i?(x a ,p a ) = i?[x a ,p a ,/i^ T (x a ,p Q ),/i^ T (x Q ,p a )], (2.6) 

where higher-order time derivatives are eliminated 
through lower-order equations of motion [llj |. corre- 
sponding to a canonical transformation. 



III. DIMENSIONAL REGULARIZATION OF 
LOCAL DIVERGENCES 

Dimensional regularization (DR) has shown its power 
in the previous 3PN calculations where unique results 



were achieved with using Dirac delta distributions as 
source functions [iJ-QjI (also see [H [lj]). Distribitional 
derivatives must still be applied but crucial for the suc- 
cess of the DR technic is the preservation of the Leibniz 
rule for differentiations and the related treatment of the 
finite part of products of singular functions. On the other 
side, at the 3PN Hamiltonian level, all occuring 1 / (<i — 3)- 
poles together with the connected logarithmic terms do 
cancel each other [12]. As efficient DR has proven in 
controlling local divergences, it failed in the treatment of 
long-range divergences. Those infrared divergences have 
nothing to do with the used point-mass model but rely 
on the very definition of the near zone. Within our for- 
malism, those divergences are still not fully under con- 
trol and need futher investigations. Only the logarithmic 
terms resulting in this context we were able to calculate 
uniquely. 

To compute the 4PN Hamiltonian i?4PN we have split- 
ted the Hamiltonian and Hamiltonian density /14PN into 
two parts, 

H iPN (d) = J d d xh APN (d) = H^ c N (d) + H^(d), 

(3.1a) 



rrloc 
-"4PN 



(d) 



d d xh l £ N (d), 



H^d)^ / d d xh}g N (d), 



(3.1b) 
(3.1c) 



where the integral over h£p N is convergent at spatial in- 
finity albeit it is locally divergent, and the integral over 
^4PN i s divergent at spatial infinity but h™p N is locally 
integrable. In the present paper we have computed the 
3-dimensional limit of the first integral. To regularize lo- 
cal divergencies we have employed the DR technic in a 
way described in detail in Ref . [12| ■ It means that in fact 
we have computed the difference 



Afi 4P N = l.imfii°p c N (d) 



d-s-3 



rrRH loc 
fl 4PN 



(3), (3.2) 



where -ff|pN loc (3) is the "local part" of the Hamilto- 
nian obtained by means of the 3-dimensional Riesz- 
implemented Hadamard (RH) regularization defined in 
Refs. The way of computing this difference was 

devised in Sees. 3 and 4 of [l2| and can be summarized 
as follows. Let us consider some global (i.e. performed 
over the whole M 3 space) integral, which develops only 
local poles. Let us denote its integrand by «(x). Then 
the value of the integral, after performing the RH regu- 
larization in 3 dimensions, usually has the structure 



7 RH (3; £l ,£ 2 ) 



A + Cl (3)(- 
+ 0(ei,e 2 ). 



hi 



ri2 
si 



c 2 (3)(l+ln^ 



(3.3) 



Here Si and s 2 are arbitrary 3-dimensional regulariza- 
tion scales. To find the DR correction to the integral 
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7 RH (3; 61,62), related with the local poles at, say, x = xi, 
it is enough to consider this part of the integrand i(x) 
which develops logarithmic singularities, i.e. which lo- 
cally behaves like 1/rf, and it is enough to consider the 
integral of this part over the ball B(x.i,£i) of radius t\ 
surrounding the particle Xi. The RH regularized value 
of this integral reads 

7^(3; £l )= / c x {?,-,n x )r^{p-Y £r x 

= c 1 (3)j\^Qy i rldr u (3.4) 

where ci(3) is the angle-averaged value of the coeffi- 
cient Ci(3;ni). The expansion of the integral I RH (3;£i) 
around S\ = equals 

jHH (3 £l) = S±W+ Cl(3 ) i n l ± + o{6 1 ). (3.5) 
£1 Si 

In the next step one computes the (i-dimensional version 
of the integral I RH (3;ei). Let us call it I\{d). It has the 
structure 

JB(xi,ei) 

= 4 {d - 3) Cl (d) [ \*- 3d rt 1 dri, (3-6) 
Jo 

where io is the scale which relates the Newtonian gravita- 
tional constant Gn with the D-dimensional gravitational 
constant Gd, 

Gd = Gn Iq 3 , 

and the number k indicates the momentum-order of the 
term [the term with k is of the order of O(p 10 ~ 2k ), where 
k = 1 , . . . , 5] . The radial integral in Eq. (|3.6p is conver- 
gent if the real part of d fulfills 5ft (g?) < 3. The expansion 
of the integral Ii(d) around e = d — 3 = reads 

h(d) = - U(3) + Cl (3) In £ + 0{6). (3.7) 

le 2 to 

Let us observe that the coefficient c' x (3) usually depends 
on lnri2, so it has the structure 

C ' 1 (3)= C ' 11 (3) + c' 12 (3)ln^. (3.8) 

Therefore the DR correction will also change the terms 
oc In rx2 ■ 

The DR correction of the integral I RH (3;£i,e 2 ) relies 
on replacing this integral by 

/ RH (3;ei,e 2 )+A/ 1 +A/ 2 , (3.9) 

where 



The corrected value of I RH (3; £1, e 2 ) thus reads 

/ R "(3; £l ,e 2 ) + A/ 1+ A/ 2 =A- Cl(3) 9 +C2(3) 

le 

-5(ci(3) + c 2 (3)) + ( C x(3) + c 2 (3)) ln^ 
= ^_ cx(3)+c 2 (3) _l (c/ii(3) + 4i(3)) 

+ (ci(3) - \c' 12 (3) + c 2 (3) - ic 22 (3)) ln^. (3.11) 

Note that all poles oc 1/ex, l/e 2 and all terms depending 
on ln^Xi In ii or lnsi, lns 2 cancel each other. The re- 
sult (|3.11[) is as if all computations were fully done in d 
dimensions. 



IV. CENTER-OF-MASS HAMILTONIAN 

Making use of the procedure described in the previous 
section we have computed DR corrections to all logarith- 
mically divergent terms contributing to H^ N (3). After 
summing up all these corrections we have obtained terms 
proportional to 1/e as well as to ln(ri 2 /4))- We have 
shown that both kind of terms can be removed from the 
Hamiltonian by adding a total time derivative. This way 
we have obtained free of poles and logarithmic terms lo- 
cal part i^p c N (3) of the 3-dimensional 4PN Hamiltonian. 
We will show here the explicit form of this part in the 
center-of-mass reference frame. Let us note that the DR 
corrections are needed only for some terms of the order 
4, 2, and in momenta. All terms of the order 10, 8, and 
6 were calculated in [6j by means of 3-dimensional RH 
regularization. 

The center-of-mass frame is defined by the condition 
Pi +P2 = 0. It is convenient to use the following reduced 
variables: r = ri 2 /(GM) (with r = |r| and n = r/r), 
p = pi/V, where M = mi + m 2 is the total mass of the 
system and /1 = mim 2 /M is its reduced mass. We also 
introduce the reduced Hamiltonian H = (H — Mc 2 )/^ 
which depends on masses only through the symmetric 
mas ratio v = /J./M (0 < v < 1/4; v = is the test-mass 
limit and v = 1/4 holds for equal masses). 

We have checked that the noncomputed part of the 
Hamiltonian, iJ"p N (3), in the center-of-mass frame and 
in the reduced form H%p N (3) = H£p N (3) / fi, is propor- 
tional to v. It contains logarithmic terms which result 
from divergences of the instantaneous near-zone metric 
when going to large distances. They were computed in 
Ref. [1] (see also Refs. @)H|). In the logarithmic terms 
s is a regularization scale [s = s/(GM), where s is a 
regularization scale with dimension of length] . 

The 4PN-accurate conservative Hamiltonian H (r, p) 
in the center-of-mass frame equals 

H{r, p) = H N (r, p) + H 1PN (r, p) + i? 2PN (r, p) 



AI a = lim I a (d) - lim Jf H (3; £1), a = 1, 2. (3.10) 

d— >3 Ei— >3 



+ J ff 3PN (r,p) +# 4 PN(r, p). 



(4.1) 
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The reduced Hamiltonians from to i?3PN can be 
found in Eq. (3.6) of [T^] (where one has to put ^static = 
and a; kinot ic = 41/24, see [HI). The 4PN reduced Hamil- 



tonian i?4PN was partially computed in Ref . @ . For con- 
venience of the reader we repeat here the main result of 
6], i.e. the Hamiltonian H4PN written in the following 
form [taken from Eq. (3.1) of [6|; here (p 2 ) = p • p and 
(np) = np]: 



c 8 H 4PN (r,p) 
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(4.2) 



The Hamiltonian (|4.2[) contains the terms C 4 i(r, p), C2i(r, p) (i = 1,2) and the coefficients C01 and C02 which were not 
computed in Ref. [6J. In the present paper we have calculated, making use of dimensional regularization, the terms 
642(1", p), C22(r,p) and the coefficient cq2- They read 



C 42 (r,p)= - 



1189789 18491 



28800 16384 



127 4035 



/ 57563 38655 



V 1920 16384 



C2»(r,p) = 



/ 672811 158177 
V 19200 ~ 49152 

1256 7403 



^ (p 2 ) + - 



21827 110099 



3840 49152 



^ (npy 



C02 



45 3072 



tt' ) (np) 4 , (4.3a) 
(4.3b) 
(4.3c) 



The still noncomputed is the coefficient C01 and the terms 
C4i(r,p) and 6*22 (r,p) of the structure 

C 4 i(r,p) = c 4 ii(p 2 ) 2 + c 4 i 2 (np) 2 p 2 + c 4 i 3 (np) 4 , (4.4a) 



C*2i(r,p) = c 2 ii(p 2 ) +c 2 i 2 (np) 2 . (4.4b) 

In Ref. @ the Poincare algebra relations were used [2(j 
to uniquely determine the three quartic in momenta and 
proportional to coefficients. In the present paper we 
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have recomputed these coefficients directly (i.e. without 
employig the Poincare algebra). 

V. CIRCULAR ORBITS 

Making use of the Hamiltonian (|4.2I) together with the 
new results of Eqs. (|4.3I) . we have computed the energy of 
binary system moving along circular orbits as a function 
of angular frequency u> (details of such computation per- 
formed at the 3PN order can be found in Ref. [19j and was 
also briefly summarized in Ref. 0). We have introduced 
the dimensionless PN parameter x = (GMw/c 3 ) 2 ' 3 and 
we have assumed (following P) that the regularization 
scale s is determined by the period P of the motion 



along circular orbits, s = cP. Making then use of the 
third Kepler's law at the Newtonian level one shows that 
s = 2ttc~ 2 x~ 3 / 2 . The 4PN-accurate binding energy of 
the system can be written in the form 



E(x; v) = — \1 + ei PN (^) x + e 2PN (^) x 2 

+ e 3PN (u) x 3 + (e 4PN (^) + ^-vlnx^j x 4 + 0(x 5 )^J , 

(5.1) 

where the fractional corrections to the Newtonian energy 
at different PN orders read 



eiPN(^) = 
e2PN(^) = 
e3PN(^) = 

e4PN(^) = 



v, 

4 12 ' 
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301 
1728 Z 



77 
31104' 



(5.2a) 
(5.2b) 
(5.2c) 
(5. 2d) 



In the above formula, the 4PN coefficient at v 2 was com- 
puted for the first time. The 4PN coefficient linear in v 
was computed numerically in Ref. @ [see also Eq. (3.1) 
in [![]. The coefficientes at v 3 and v 4 computed for the 
first time in Q was recently confirmed by independent 
calculation made in Q. 

An important feature of the two-point-mass dynam- 
ics along circular orbits is the existence of the last sta- 
ble circular orbit (LSCO) [2l|. In the test- mass limit 
v = 0, the LSCO occurs for xlsco = 1/6 what cor- 
responds to the minimum of the function E{x;v = 0). 
Therefore the most straightforward way of locating the 
LSCO for v > relies on looking for minimum of the 
function E(x; v). In Fig. Q] (left panel) we have depicted 
the location £lsco of the LSCO as function of the sym- 
metric mass ratio v for successive PN approximations 
from 1PN up to 4PN. For a given value of v, the loca- 
tion of the LSCO was obtained by numerically solving 
the equation dE(x; u)/dx = 0. Let us note that the 4PN 
prediction for the location of the LSCO in the test-mass 
limit is 0.179467, ^7.7% larger than the exact result. 
We have obtained the following locations of the LSCO 
in the equal- mass case (y = 1/4) for the approximations 
from 1PN up to 4PN: 0.648649 (1PN), 0.265832 (2PN), 
0.254954 (3PN), 0.236597 (4PN). 

We have also found the location of the LSCO by means 
of more refined method, namely the "j-method" intro- 



duced and used at the 3PN level in Ref. [2l[. In this 
method the xlsco is defined as the minimum of the func- 
tion j 2 (x), where j = J jiGmim-i) is the reduced angular 
momentum of the system (JJ is the total angular momen- 
tum). Moreover, motivated by the form of the function 
j 2 {x) in the test-mass limit, j 2 (x) = l/(x(l — 3a:)), Pade 
approximants are used instead of direct Taylor expan- 
sions (all used approximants have a pole for some x po i e (u) 
which is related with the test-mass "light-ring" orbit oc- 
curing for x\ x = 1/3 in the sense that x po \ e {v) — s- 1/3 
when v 0). The 4PN-accurate function j 2 (x) has the 
symbolic structure (l/x)(l + x + . . . + x 4 + x 4 \nx). In 
the j-method the Taylor expansion at the 1PN level of 
the symbolic form 1 + x is replaced by Pade approximant 
of type (0,1), at the 2PN level 1 + x + x 2 is replaced by 
(1,1) approximant, at the 3PN level 1 + x + x 2 + x 3 is re- 
placed by (2,1) approximant, and finally at the 4PN level 
1 + x + x 2 +x 3 +x 4 is replaced by (3,1) Pade approximant 
[the explicit form of the (0,1), (1,1), and (2,1) approxi- 
mants can be found in Eqs. (4.16) of [HI]. The results 
are illustrated in Fig. [1] (right panel). At all PN levels 
the test-mass result is recovered exactly. We see that the 
curves corresponding to 3PN-accurate and 4PN-accurate 
calculations almost coincide. The locations of the LSCO 
in the equal- mass case v = 1/4 for the approximations 
from 1PN up to 4PN are as follows: 0.162162 (1PN), 
0.185351 (2PN), 0.244276 (3PN), 0.247515 (4PN). 
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FIG. 1. The location £lsco of the LSCO as a function of the symmetric mass ratio v for successive PN approximations from 
1PN up to 4PN. Solid curves correspond to 1PN, dashed curves to 2PN, dotted curves to 3PN, and dashed/dotted curves to 
4PN. Left panel: 2;lsco is computed as the minimum of the function E(x; v). Right panel: xlsco is obtained by means of the 
j-method. 
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